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ABSTRACT 

We present an updated cumulative size distribution (CSD) for Jupiter Family comet 
(JFC) nuclei, including a rigourous assessment of the uncertainty on the slope of the 
CSD. The CSD is expressed as a power law, 7V(> rjsr) oc r N q , where tn is the radius 
of the nuclei and q is the slope. We include a large number of optical observations 
publis hed by ourselves a nd others since the comprehensive review in the Comets II 
book (jLamv et al.ll2004l ). and make use of an improved fitting method. We assess 
the uncertainty on the CSD due to all of the unknowns and uncertainties involved 
(photometric uncertainty, assumed phase function, albedo and shape of the nucleus) by 
means of Monte Carlo simulations. In order to do this we also briefly review the current 
measurements of these parameters for JFCs. Our final CSD has a slope q = 1.92±0.20 
for nuclei with radius tn ^ 1.25 km. 

Key words: comets. 



1 INTRODUCTION 

The size distribution of any population of solar system small 
bodies is of critical importance in constrai ning their forma- 
tion and subsequent collisional evolution. iDohnanvil (|l969l ) 
showed that a collisionally relaxed population of self-similar 
bodies, with the same strength per unit mass, has a charac- 
teristic power law size distribution with a slope of 2.5. For 
a colli sional population of gravity controlled (strengthless) 
bodies lO'Brien fc Greenberel (l2003h demonstrated that the 
expected size distribution has a shallower slope, 2.04. Hap- 
pily, the size distribution is also one of the more straight- 
forward characteristics of the population to determine, as at 
least reasonable estimates of the size of bodies can be made 
with snap-shot observations. Time-series data allow bet- 
ter measurements of their sizes, as this removes uncertain- 
ties due to their rotational light curve. For Jupiter Family 
comets (JFCs) these observations are generally made when 
the comet is at a large distance from the Sun, and there- 
fore more likely to be inactive, so a brightness measurement 
for the bare nucleus can be made. Converting the measured 
optical magnitude to the size of an object depends on its 
albedo and phase function. There are only 7 JFCs for which 
both of these are independently well measured (2P/Encke, 
9P/Tempel 1, lOP/Tempel 2, 19P/Borrelly, 28P/Neujmin 1, 
67P/Churyumov-Gerasimenko and 81P/Wild 2; three of 
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which have a size measurement from resolved imaging by 
spacecraft [see Tables [3] and U for references] ) . A few more 
comets have measurements of one or the other. Where these 
parameters are not known, we are forced to assume values 
for them; typically 4% for the albedo and a linear phase func- 
tion of P = 0.035 mag deg -1 . Previous measurements of the 
size distribution of JFCs have followed similar assumptions. 

The distribution is generally plotted (fig. [1} as a cumu- 
lative size distribution (CSD), expressed in terms 



N(> m) oc rj 



(1) 



N is the number of nuclei with radius, rn : larger than 
FN, and q is the slope. A number of measurements of 
these distribution co-efficients have been made, generally 
based on snap-shot observations of a large number of nuclei. 
lLowrv et ail (120031 ) est imated q = 1.6 ±0.1 based on a sam- 
ple of 33 comets, and IWeissman fc Lowrvl (|2003t ) updated 
this to q = 1.59 ± 0.03, based on 41 JFCs with r N ^ 1.4 
km. They chose to fit the size distribution to only those 
nuclei with tn 1? 1.4 km as the slope of the CSD is ap- 
proximately constant above this radius, while below it there 
is a sharp cut-off. This break may imply a relative paucity 
of small nuclei compared with the expected n umber from a 
contin uati on of the power law f rom l arger sizes; iMeech et al.l 
(|2004h and lFitzsimmons et al.l (120111 ) show that such a break 
is most likel y real and not due to observational biases (as 
suggested by lLamv et al.l (|2004l )) by modelling realistic ob- 
servational surveys. A more recent size distribution from 
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Table 1. Previous JFC size distribution estimates. 



Reference 



AW-Ntot r N range 



Fern andez et al, (1999 ) 
Lowrv et al. (20031 
Wcisjmari_^IjOwrv i (2003) 
Meech et al. (2004) 
Meech et al. (2004) 



Lamv et al. (2004) 
Weissman et al. (2009) 
Tancredi et al. (2006) 



12/64 
16/19 
41/54 
38/48 
21/48 
29/65 
41/67 
32/72 



2.1 - 
1.4 - 
1.4 - 
1.0 - 
2.0 
1.6 - 
1.4 - 
1.7- 



3.3 
3.6 
15 
10 
- 5 
15 
6.0 
4.5 



2.7 ±0.3 

1.6 ±0.1 
1.59 ±0.03 
1.45 ±0.05 
1.91 ±0.06 

1.9 ±0.3 
1.94 ±0.07 

2.7 ±0.3 



IWeissman et al.1 (120091 ) has a value of q = 1.94 ± 0.07 based 
on 41 JFCs with tn> 1.4 km. Another estimate comes from 
iFernandez et al.l (j 19991 ). who used selected data in quality 
classes 1-3 (uncertainties on my (1,1,0) up to ±1 mag.) 
from the catalogue presented bv lTancredi et all (|2000l ). with 
cut-offs in both absolute mag nitude and perih elion dis- 
tance. The discrepancy between Fernandez et alj 's estimate 
of q — 2.65 ±0.25 and those of lLowrv et al.l can be explained 
by these cut-offs (leaving only 12 comets) an d the large 
uncertainties on mag nitudes in iTancredi et alj 's catalogue. 
I Tancredi et al.l (|2006h presented a n updated ca t alogu e and 
find q = 2.7 ± 0.3 for r N 1.5 km. iMeech et all (|2004l ) esti- 
mate q = 1.45 ± 0.05 over the range 1 ^ rN ^ 10 km, and a 
steeper q — 1.91 ±0.06 in the range 2 ^ tn ^ 5 km, showing 
the la rge dependence on the choice of size range. iHicks et ahl 
|2007l ) estimate q = 1.50 ± 0.08 from Near Earth Asteroid 
Tracking (NEAT) survey data, although this makes use of 
observations of active comets and a coma subtraction tech- 
nique, so is not included with th e other surveys of d istant in- 
active comets in TableED Finallv. lLamv et al.l (|2004T > collated 
the data from most of these catalogu es, together with their 
own unpublished results, those from iLicandro et al.l (|2000l ) 
and also from other papers on individual comets, to calcu- 
late q — 1.9 ± 0.3 for JFC nuclei with tn ^ 1.6 km. 

These previous estimates are listed in Table [TJ where 
we list the number of comets included in the fit (and the 
total number considered in the survey), the range in sizes 
over which the authors fit the linear part of the CSD, and 
the resulting slope q. The uncertainty on q is that quoted 
by each author, and is generally the formal uncertainty from 
a least squares (or similar) fit to the line, despite the fact 
that technically one cannot use such a fitting technique as 
the points in a cumulative distribution are not independent. 
These previous works also make no attempt to assess the 
uncertainty on q contributed by the assumptions on phase 
function and albedo or the uncertainty from the photometry. 
In this paper we present an updated size distribution based 
on new photometry of distant JFC nuclei, using a censored 
data analysis technique to produce the CSD, and make a 
rigourous assessment of the uncertainty on q due to all the 
unknown factors. 



2 METHOD 

2.1 Censored data analysis 

The majority of previous measurements of the comet nu- 
cleus CSD follow a similar procedure. For each nucleus a 
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Figure 1. Reference CSD with our data set and the usual as- 
sumptions. This shows the normal shape of the CSD, with a lin- 
ear part and a cut-off. The error bars shown are calculated using 
the Kaplan-Meier statistic. Each fit within a MC run produces a 
CSD like this with a slightly different distribution of points. 



best estimate of the intrinsic radius is arrived at either from 
the investigators own data or though consideration of sev- 
eral reported detections. The resulting radii are combined to 
form a CSD and a straight line is fit through the CSD using 
linear regression. Although this is relatively straightforward, 
there are intrinsic problems to this method. Aside from the 
incorrect use of slope uncertainties from linear regression 
techniques mentioned above, a more important drawback 
is that it ignores non-detections or upper limits to nuclear 
radii. Not only is this is an important source of information, 
ignoring these non-detections can act to instill significant 
bias. For example, consider a single imaging survey of inac- 
tive nuclei with a well defined limiting magnitude for detec- 
tion. Nearby comets might all be detected, but the fraction 
of comets observed would decrease as their geocentric and 
heliocentric distances increased. The resulting CSD of nu- 
clei would contain contain a progressively smaller fraction 
of the true population at smaller sizes, giving an observed 
CSD with a shallower power-law distribution than in real- 
ity. In the past, breaks in observed luminosity distributions 
have often been interpreted as due to the onset of such in- 
completeness. To our knowledge, the only pr evious study 
that h as accounted for upper limits was that o f lMeech et al.l 
|2004h . 

Many fields of astronomical research face a similar prob- 
lem of censored data, where a sample contains both di- 
rectly measured values and upper limits. The combina- 
tion of ob s erved and censored data was first tackled by 
lAvni et al.1 (|l980l ) via a maximum likelihood method. Estab- 
lished statistical methods as applied to cens ored astr o nomi - 
cal d atasets were subsequ e ntly d escribed bv lSchmittl (|l985l ) 
and iFeigelson fc Nelson! i|l985h . These authors showed 
how the cumulative distribution function of a censored 
dataset can be constrai ned via the Kaplan-Meier estimator 
|Kaplan fc Meier] 1 19581 ). This statistic normally estimates 
the probability of measuring a value for an observable less 
than or equal to a given value at the non-censored values of 
the dataset, or in our case for radii of R km, P(rjv ^ R). 
As long as the observed sample is randomly picked from the 
population, then this will be equivalent to estimating the 
relative number of comets with radii ^ R. Inverting this 
function then gives the normal form of the CSD. 
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For cometary nuclei at large heliocentric distances the 
probability of non-detection will be larger as nucleus size 
decreases, hence at first sight the sample will be biased to- 
wards large nuclei. However, other effects act to randomise 
the censoring. First, comets are observed at various positions 
in their orbits, which when combined with different sized 
telescopes, imaging cameras and exposure times leads to a 
wide variation in effective detection limits. As an example, a 
600-second exposure at R-band with the 4.2m William Her- 
schel Telescope can detect a point source at mu = 24.5 at a 
signal-to-noise of ~ 6. Assuming an optical albedo of 0.04, 
this allows direct detection of a rjv = 0.8km nucleus at op- 
position at Rh = 5.2 AU, or the same nucleus at opposition 
at A = Rh = 4AU but at a phase angle of ~ 20°. Second, 
when comets become active the shielding effect of their dust 
comae can result in censoring of even large cometary nuclei. 
As the strength and onset of activity appears to strongly 
vary from comet to comet, this further randomises whether 
an inert nucleus will be directly imaged or not. In this study 
we therefore assume that the observed cometary nuclei ap- 
proximates a randomised sample of the JFC population. 

The use of the Kaplan-Meier statistic is accurate as long 
as the quantity being measured and the censoring variable 
(here the nuclear radius and the limiting magnitude) are in- 
dependent. The limiting magnitude is set by the instrumen- 
tal setup, exposure times, and activity level of the comet, 
while the individual comets surveyed will normally depend 
on their orbital visibility at a particular date. Therefore this 
requirement should be true of most cometary observations 
where non-detection is due to large distance leading to a 
faint apparent magnitude. For nucleus measurements cen- 
sored by activity, it will also hold. Although weak corre- 
lations of active area with orbital parame ters have previ- 
ously been suspected (|A'Hearn et al.|[l995l ). t here appears 
to be no link between nuclear radius and orbit (|Lamv et al.l 
l2004h . and the orbital positions of comets are randomised 
by the dates of observation. Finally, an important aspect 
of the Kaplan-Meier statistic is that in principle it allows 
an estimate of the uncertainty of the luminosity function at 
each value where an uncensored measurement is made. One 
method of calculating this is giv en in equations (11) and 
(12) of lFeigel son fc Nelson! (|l985h . which we adopt in this 
work. We therefore plot the CSD in terms of P(rN > R) 
with uncertainties at each point from this statistic, which is 
equivale nt to N(r^ > R) and give s the slope q we are inter- 
ested in. iFitzsimmons et all l|201lf ) give more details on this 
method and its application to comet size distributions. 



2.2 Fitting the CSD 

We assume that the CSD is well described by two power laws 
(i.e. straight lines on a log-log plot), as we expect a break in 
the power law at small sizes. The power law describing the 
CSD at larger sizes is the one of primary interest, and is the 
one we refer to when discussing the slope q elsewhere in this 
paper, as the data at smaller sizes suffers more from obser- 
vational incompleteness. The choice of radius to define this 
break is important, as some of the variation of slopes found 
by the various papers listed in Table [1] can be explained by 
the different radius ranges over which the authors chose to 
fit the CSD. We attempt to independently choose the cor- 
rect break radius by testing all possible radii in the range 
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Figure 2. Combined \ 2 f° r two power law fits either side of a 
cut-off, as a function of the selected cut-off radius, for the refer- 
ence CSD shown in fig. [T] We select the cut-off radius with the 
minimum \ 2 as the best description of the overall shape of the 
CSD, and therefore the corresponding q value from the fit to the 
CSD at larger radii as the best slope for that particular CSD. In 
this case, the minimum is at a radius of 1.26 km. 



of the comet sizes. We fit two straight lines on either side 
of the break, and compute the \ 2 statistic for the fit of this 
model to the CSD, for each test radius (fig. [2)1. The mini- 
mum X 2 ( r ) then gives us the optimum break radius, which 
tends to fall around 1.2 km. 

While strictly speaking a \ 2 fit is not appropriate in this 
case, as a cumulative distribution does not have independent 
points, our goal is to find a simple description of the fit 
that gives a reasonable measure in an individual case. The 
formal uncertainty on each fit is not used, since we find 
the overall uncertainty from the standard deviation around 
the average of all fits. We investigated the results from a 
number of fitting techniques, to be sure that the result from 
the x 2 technique is representative. We performed these tests 
on our reference size distribution (fig. [1] see section [3| ■ The 
X 2 minimisation technique gives a value of q = 1.97 for this 
data set. 

One common approach to fitting a line to a cumulative 
distribution is to use a maximum-likeli hood fit. Using the 
formulae given by iBadenes et al.l (|2010l) for the maximum 
likelihood power law index, we find q — 1.95 ±0.29, a nearly 
identical result to the x 2 method. 

An alternative, as used bv lLamv et alJ (|2004l ) . is to test 
the observed distribution against various model distribu- 
tions using a Kolmogorov-Smirnov test. In this case the best 
fit is the one for which the probability that the samples are 
drawn from the same distribution is the highest. Applying 
this approach to our test size distribution returns a best 
fit value of q = 2.03 with a probability Pks of 0.99, for the 
same cut off radius. This method has a significant draw back 
though, as the Kolmogorov-Smirnov test is not valid when 
the parameters of the model (in this case the cut-off radius 
and normalisation of the power law) are estimated from the 
data, which is necessary in this work. 

A further disadvantage with both the maximum- 
likelihood and K-S techniques is that they fail to take advan- 
tage of the fact that the K-M statistics return uncertainties 
on each radius, which include information from upper limits 
and should be taken into account in an ideal fit. On the other 
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hand, while the fit is not appropriate due to the fact the 
points are not independent, the use of K-M statistics mean 
that the error bars give proper weighting to each point. In 
conclusion, while none of these statistical methods are fully 
appropriate to these data, the fact that they all result in a 
similar value for q gives us confidence that the result from 
X 

the size distribution over this size range 



2 is a reasonable estimate of the power-law exponent for 



2.3 Monte Carlo techniques 

The effective radius of the nucleus is given by 

A R r 2 N = 2.238 x 10 22 i? 2 A 2 10 - 4(mo ~ m « +,3a) (2) 

where Ar is the geometric albedo, R^, A and a are the 
heliocentric and geocentric distances (AU) and phase angle 
(degrees) at the time of observation, j3 is the linear phase 
coefficient (mag deg -1 ) and tur and mg = —27.09 are the 
apparent m agnitudes of the comet and the Sun, both in 
the i?-band (|Russelllll916h . The orbital position parameters 
(Jih, A, a) for a given observation and the magnitude of 
the Sun are known to a high level of precision, but the re- 
maining parameters all contribute to the uncertainty on the 
reported radius. In addition, this equation converts a cross- 
sectional area (found as the area is directly proportional to 
the reflected flux) to an effective radius assuming a spheri- 
cal nucleus, while the shape is also generally unknown (see 
section |4~4]) . We examine the effect of each of these sources 
of uncertainty by applying a Monte Carlo (MC) technique. 
For each comet observation in our data set (Table [2]| we se- 
lect a random value for each of the unknown parameters (/?, 
An) from a distribution of possible values and a magnitude 
within the range turgor, and use this to generate a radius. 
We correct this radius for the effects of the non-spherical 
shape of the nucleus by selecting an ellipsoidal shape from 
a suitable distribution of axial ratios. We then generate a 
CSD by fitting the resulting radii as described in the previ- 
ous subsections. We repeat this procedure many times (gen- 
erally, N = 1000), selecting a different random value for 
the parameters for each comet each time, which allows us 
to measure the effect of varying the parameters within our 
chosen distributions. We then measure the average value of 
< q > found from the N different CSDs; the standard devi- 
ation on this a(q) then incorporates the uncertainty on the 
varied parameters. In section [4] we vary each parameter in 
turn (while keeping the other unknown values fixed at de- 
fault values) to assess the relative contribution to the total 
uncertainty of each parameter. Finally we allow all of the 
parameters to vary simultaneously to measure the overall 
uncertainty on q. 



3 DATA SET 

The data used for this investigation are selected from all 
JFC nucleus photometry found in the literature (only prop- 
erly calibrated professional data sets), including our own 
previously published work. The majority come from the sur- 
veys listed in the introduction that made their own esti- 
m ates of the CSD, w hich were all collated in the review 
bv lLamv et al.l (|2004l ). We also include the many observa- 
tions published in the six years since the publication of that 



work, which have added new comets to the database and im- 
proved the accuracy of the radius measurements of others 
(in many cases measuring light curves where previously only 
snap-shot measurements had been taken ). We include the 
surveys of Lowrv fc Fitzsimmonsl (l2005f). Snodgrass et al.l 
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ey s ot IL/Owrv &: i 1 ltz simmons llzltuoll. la nodgrass et al.l 
a 120061. l200Sf). iMazzotta Epifani et al] (|2007t I2008T ) 



and Iweissman et al] (120081). and also papers on i ndivid- 
ual co mets from iFernandez et al.l (l2005h [2P/Encke1 , Ijewittl 



(|2006T ) [D/1819 Wl (Blanpain)l ISnodgrass et al.l (|2010h 
[103P/Hartley 2], iTubiana et"afl (120081) [ 67P/Churvumov- 
Gerasimenko] and Pittichova et al] (|2008l ) [21P/Giacobini- 



Zinner]. The references for each comet are given in Table 

m 

To apply our MC method, it is necessary that we choose 
one observation at a given Rh, A and a for each comet, and 
use equation [2] to generate the corresponding radius for a 
chosen set of parameters. For some comets there have been 
multiple observations at different epochs, and we therefore 
require suitable criteria to select the 'best' observation. We 
must do this rather than attempting to take an average value 
of the radiu s for the comet fr om separate observations (as 
was done bv lLamv et al.l (|2004l ) , for example) as such a com- 
bination would require an assumed phase function to correct 
for the different orbital positions, and the effect of the un- 
known phase function is one of the parameters we investi- 
gate. 

The criteria we use are: 

(i) Is there a detection of the inactive nucleus? If not, we 
use the strongest upper limit, and flag the value as a limit 
(except in the cases where space based observations allowed 
an accurate coma subtraction). 

(ii) Is there a light curve? Where there is a direct nucleus 
detection we prefer light curves over snap-shots, as this re- 
moves some of the uncertainty due to the unknown shape. 

(iii) Finally, where we are left with a choice between snap- 
shot observations, we use the result with the smallest pho- 
tometric uncertainty and smallest phase angle to minimise 
uncertainty. 

The final list of selected observations is given in Table [2] 
The table also includes the radii calculated under the stan- 
dard assumptions (/3 = 0.035 mag deg -1 and Ar = 0.04 for 
all comets; photometric values are precise and correspond 
directly to the size of a spherical nucleus). This is for ref- 
erence only as the radius found for each comet varies in 
each MC run, over a range set by the variable parameters, 
throughout the rest of the paper. The data set contains 86 
comets. 71 are nucleus detections, 15 are upper limits from 
either active or non-detected comets. 21 have light curves, 
and another 2 have a constraint on a magnitude range but 
not a full light curve. 12 have known j3 (although these are 
often very rough estimates), 11 have known albedo and 7 
have both of these parameters measured. 

For reference, fitting a CSD to this data set without 
varying any parameters (i.e. following the procedure ap- 
plied to previous works; taking the reference radii based on 
the standard assumptions listed in Table [2|) gives q — 1.97, 
with a cut off at 1.26 km (fig.[l|. For comparison, using the 
measured phase functions and albedos for the comets where 
these values are known (and using the assumed values given 
above for other comets) gives q = 2.00, and the same cut off. 
There are relatively few comets with known j3 and/or Ar 
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Table 2. Photometry database used for input 
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Reference 



Fernandez ct al. (2005) 
Lamv et al. (2009) 
j Lowr^^Weissman (2003) 
Snoderass et al. (2005) 
Lamv et al. (2007) 
Jewitt fc Luu (1989) 
Snoderass et al. (2005) 
Snoderass et al. (2006) 
Lamv et al. (1998) 
Pittichova et al. (2008) 
^owr^^^Veissman (2003) 
Boehnhardt et al. (1999 ) 
Delahodde et al. (2001) 




Meech et al. (1993) 
Luu fc Jewitt (1992) 
Snoderass et al. (2008) 
Lamv et al. (2009) 
Snoderass et al. (2008) 
^owrv_^J^itzsimmons 
Snoderass et al. (2008) 
Lowrv et al. (2003) 
Boehnhardt et al. (2002) 
Snoderass et al. (2006) 
Jewitt fc Sheppard (2004) 
Lowrv et al. (2003) 
Lamv et al. (2009) 
JjOwrv_fc_J^itzsimmons 
Licandro et al. (2000 ) 
Meech et al. (20041 
^owrv_fc_^itzsimnioris 
Meech et al. (2004) 
^owrv_^J?itzsimmons (2001) 
Lamv et al. (2009) 
Lowrv et al. (2003) 
Lamv et al. (2009) 
Licandro et al. (2000 ) 
^owrv_^_J?itzsimmons (2001) 
Tubiana et al. (2008) 
Lowrv et al. (1999) 
Snoderass et al. (2008) 
Lamv et al. (2009) 
Boehnhardt et al. (1999) 
Lamv et al. (2004) 
Weissman et al. (2008) 
Lamv et al. (2004) 
I^wr^t^eissman (2003) 
Lowrv et al. (1999) 
Weissman et al. (2008) 
Lamv et al. (2004) 
Lamv et al. (2009) 
Meech et al. (2004) 
Lamv et al. (2004) 
Snoderass et al. (2005) 
Snoderass et al. (2008) 
Lowrv et al. (2003) 
Mazzotta Epifani et al. (2008 ) 



a The observed variation in magnitude, corresponding to a minimum axial ratio. 

^ The radius (km) calculated for a fixed j3 = 0.035 mag deg - 1 and albedo of 4%, for reference. 

c We include 2P/Encke, which is an ecliptic comet although not technically a JFC by the Tisserand parameter definition. 
d From reported ^-magnitude and (V - R) = 0.50 ± 0.01 from lLi et all d2007al) 
e From reported V-maenitudc and assumed (V — R) = 0.50. 
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Table 2. Photometry database used for input (continued) 



Comet 



Rcfcrencc 



lOOP/Hartley 1 4.84 4.06 

103P/Hartley 2 5.60 4.70 

104P/Kowal 2 3.94 3.31 

106P/Schuster 1.67 0.78 

HOP/Hartley 3 3.60 2.66 

lllP/Helin-Roman-Crockett 3.47 2.53 

112P/Urata-Niijima 2.30 1.50 

114P/Wiscman-Skiff 3.75 2.95 

115P/Maury 6.38 5.43 

116P/Wild 4 4.72 3.82 

117P/Hclin-Roman-Alu 1 3.29 2.56 

118P/Shoemaker-Levy 4 4.71 3.72 

120P/Mueller 1 3.89 3.05 

121P/Shoemaker-Holt 2 3.92 3.43 

128P/Shoemaker-Holt 1 4.99 4.16 

131P/Mueller 2 3.48 2.65 

136P/Mueller 3 4.83 4.04 

137P/Shoemaker-Levy 2 6.95 6.17 

143P/Kowal-Mrkos 3.40 2.47 

147P/Kushida-Muramatsu 2.83 2.30 

160P/LINEAR 3.98 3.41 

179P/Jedicke 5.52 4.71 

243P/NEAT 3.97 3.68 

D/1819 Wl (Blanpain) 1.64 0.72 

P/2001 H5 (NEAT) 4.61 3.62 

P/2002 T5 (LINEAR) 5.24 4.46 

P/2003 SI (NEAT) 3.54 3.05 

P/2004 D029 (Spacewatch-LINEAR) 4.22 3.33 

P/2004 H2 (Larsen) 3.71 3.04 

P/2004 H3 (Larsen) 3.71 2.97 



8.20 22.06 0.06 2.48 Weissman et al, (2008) 

4.50 ^24.50 0.04 ^ 1.02 Snodgrass et al. (20101 

12.66 23.05 0.93 1.12 Lowrv et al. (20031 

23.00 18.91 0.03 0.89 Lamv et al. (2009) 

6.58 20.58 0.10 2.33 Weissman et al. (20081 

5.60 21.48 0.05 1.39 Mazzotta Epifani et al. (20071 

19.20 20.96 0.04 0.86 Lamv et al. (20091 

10.90 22.95 0.11 0.97 Snodgrass et al. (20081 

2.56 24.86 0.07 1.10 Meech et al. (20041 

6.00 ^21.36 0.06 ^3.03 Weissman et al. (20081 

13.50 ^16.31 0.03 < 16.36 Mazzotta Epifani et al. (20081 

3.50 21.54 0.20 2.61 Lowrv et al. (20031 

8.80 23.52 0.13 0.77 Snodgrass et al. (20081 

13.50 20.77 0.01 0.15 3.34 Snodgrass et al. (20081 

7.29 22.02 0.09 2.63 Lowrv fc Weissman (20031 
10.10 22.77 0.09 0.87 Snodgrass et al. (20081 

8.00 23.70 0.20 1.15 Mazzotta Epifani et al. (20081 

5.30 22.86 0.03 3.57 Snodgrass et al. (20061 
8.20 18.52 0.02 5.42 Jewitt (20021 

18.70 25.48 0.10 0.46 0.20 Lamv et al. (20041 

12.50 23.69 0.17 0.87 Snodgrass et al. (20081 

6.90 22.63 0.07 2.47 Snodgrass et al. (20081 

14.30 22.70 0.30 1.52 Mazzotta Epifani et al. (20081 

20.70 22.40 0.10 0.16 Jewitt (20061 

3.13 23.29 0.11 1.10 Lowrv fc Fitzsimmons (20051 

7.50 >19.10 0.30 sg 11.40 Mazzotta Epifani et al. (20081 

15.20 ^21.37 0.03 ^2.10 Mazzotta Epifani et al. (20081 

7.00 >20.13 0.03 ^ 4.23 Mazzotta Epifani et al. (20081 

13.10 ^21.59 0.01 < 1.91 Snodgrass et al. (20081 

12.00 24.19 0.21 0.55 Snodgrass et al. (20081 



(see Tables [3] and [4]), and these tend to be the larger comets 
which have a weak affect on the CSD, so the difference is 
not large. 

4 VARIATION OF PARAMETERS 

4.1 Photometric uncertainty 

The first source of uncertainty in the radius, and the most 
straight forward to account for, is the uncertainty an in 
the measured magnitude of the comet uir. This includes 
the Poisson noise from the CCD photometry, uncertainty 
from the standard star calibration, etc. The true value can 
fall anywhere within the Gaussian distribution centred on 
the reported value, with the standard deviation set by the 
reported error bar. In most cases an ^ 0.1 mag, and this is 
not the most significant uncertainty on tn- The result of a 
MC run allowing the comet magnitudes to vary within the 
Gaussian distributions defined by the reported mn, an for 
each, and holding all other parameters fixed at the default 
values given above, gives q = 1.93 ± 0.10 and an average 
cut-off of 1.2 km. 

4.2 Phase function 

For all previous size distributions, a constant Solar phase 
function has been assumed for all comets. The work pre- 



Table 3. Phase function measurements for JFCs. 



Comet /3 Reference. 



2P 


0.060 ± 0.005 


Fernandez et al. (20001 




0.049 ± 0.004 


Boehnhardt et al. (20081 




0.053 ± 0.003 


Weighted mean 


9P 


0.046 ± 0.007 


Li et al. (2007a1 


10P 


0.037 ± 0.004 


Sekanina (19911 


19P 


0.043 ± 0.009 


Li et al. (2007b1 


28P 


0.025 ± 0.006 


Delahodde et al. (20011 


36P 


0.060 ± 0.019 


Snodgrass et al. (20081 


45P 


~ 0.06 


Lamv et al. (20041 


47P 


0.083 ± 0.006 


Snodgrass et al. (20081 


48P 


0.059 ± 0.002 


Jewitt & Sheppard (20041 


67P 


0.076 ± 0.003 


Tubiana et al. (20081 


81P 


0.0513 ± 0.0002 


Li et al. (20091 


143P 


0.043 ± 0.001 


Jewitt et al. (20031 



sented in this paper was motivated in part by a simple 
demonstration that changing the assumed value changed the 
resulting value o f q: increasing the assumed j3 increases q 
(|Snodgrassll2006j ). To investigate the more realistic situa- 
tion of P varying between comets within some distribution 
we developed the MC technique presented here. We recreate 
the earlier result when choosing the simplest distribution; a 
constant value for all comets. 

To choose an appropriate distribution we briefly re- 
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Table 4. Albedo measurements for JFCs. 







0.02 



0.08 



0.04 0.06 
/? (mag deg ) 
Figure 3. Distribution of measured phase functions for JFCs. 



view the measured values. Phase functions for comets are 
difficult to measure, since it is necessary to obtain light 
curves at multiple epochs to remove rotational variations 
from the photometry and leave the variation due to chang- 
ing phase angle. Even relatively well measured phase cu rves, 
such as that for 28P/Neujmin 1 jDelahodde et al.ll200lh , are 
still quite uncertain, and rough estimates from only a few 
epochs are v ery approximate even if formally good fits to 
the data (e.g. ISnodgrass et al.ll2008l '). The best results come 
from spacecraft data, although we note that spacecraft may 
measure phase functions to high a, while ground based mea- 
surements generally only cover a < 15°. We list all JFCs 
with estimates in Table[3] and plot the distribution of values 
in fig. It is noticeable that the measured phase functions 
are generally steeper than the normally assumed values of 
(3 = 0.035 or 0.04 mag deg -1 ; the mean value is 0.053 ± 
0.016. This is due to a number of recent measurements of 
steep phase functions for JFCs, while the larger (brighter, 
easier to observe) nuclei observed previously all had shal- 
lower slopes. We tested for any correlation between size of 
nucleus and phase function, but none exists. The j3 values 
are approximately normally distributed around the mean, 
so we use a Gaussian distribution for input into the MC 
simulations. 

A MC run with a distribution matching this observed 
one gives q = 1.95 ± 0.07. Increasing the mean /3 increases 
q, and decreasing it decreases q, as expected. Changing the 
width of the distribution also affects the CSD slope, a nar- 
rower distribution (/3 = 0.053 ± 0.005) gives a steeper slope 
and a smaller uncertainty, q — 1.99 ± 0.02, while a wider 
distribution (a — 0.025) gives a shallower slope and larger 
uncertainty (q = 1.91 ± 0.11). 



4.3 Albedo 

The conversion from a reflected flux to a radius requires 
knowledge of the albedo, which generally must be assumed 
for comets. Under the normal assumption that all comets 
have the same albedo the CSD slope is directly proportional 
to the slo pe of the 'luminosi ty' (reflected flux) distribution 
(?s = 5ar.. llrwin et al.l (1 19951 )). and so a change of assumed 
albedo will not affect q. However, this does not hold for a 



0.1 



Comet 



Reference 



21' 

9P 



10P 
19P 



22P 

28P 

49P 

67P 

81P 

103P 

162P 



0.050 ± 
0.064 ± 
0.046 ± 
0.072 ± 
0.061 ± 
0.030 ± 
0.029 ± 
0.072 ± 
0.033 ± 
0.048 ± 
0.03 ± 
0.045 ± 
0.054 ± 
0.064 ± 
0.028 ± 
0.037 ± 



0.030 
0.013 
0.015 
0.016 
0.008 
0.012 
0.006 
0.020 
0.006 
0.010 
0.01 
0.019 
0.006 
0.010 
0.009 a 
0.014 



Fernandez et al. (2000) 
Li et al. f 2007a) 
Lisse et al. (2005 ) 
Fernandez et al. (2003) 
Weighted mean 
A'Hearn et al. (1989) 
Buratti et al. (2004) 



Weighted mean 
Lamv et al. (2002) 
Jewitt fc Meech (1988) 
^jamrjms_et_aL (1995) 
Kellev et al. (2009) 
Li et al. (2009) 
Lisse et al. (2009) 
Fernandez et al. (2006) 



a Uncertain due to di fficulties with the optical photometry, see 
ISnodgrass et al.l i20ld ) 



0.02 0.04 0.06 

Albedo (R-band) 
Figure 4. Distribution of observed i?-band albedos. 



distribution of albedo values that vary between comets, so 
we use our MC technique to assess the affect on the CSD of 
a variable albedo. 

The distribution of albedos for JFC nuclei is poorly con- 
strained. Only when the radius can be determined indepen- 
dently of the flux can the albedo be found, and this has been 
done for very few comets (generally by simultaneous obser- 
vation in the optical and thermal infrared, where a thermal 
model of the nucleus returns the size for comparison with the 
reflected optical flux, or from spacecraft observations). We 
list all published JFC albedos in Table [4] converting albe- 
dos to .R-band where necessary, and plot the distribution of 
_R- band values in fi g. 3] W e note that the values obtained 
bv lLi et all (|2007d lb[ [2009l '). from modelling of disk resolved 
photometry from spacecraft data, are higher than all but 
one of the other measurements. This suggests some under- 
lying problem in comparing results from this technique with 
radiometric methods, but we take them at face value for 
the purposes of this study. The mean value is 0.044 ± 0.013, 
close to the usually assumed 4%. The small number of known 
values means that the shape of the distribution if not well 
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defined; we choose to use a Gaussian distribution as this 
seems reasonable. 

For a Gaussian with the observed mean and standard 
deviation, the CSD has q — 1.85 ± 0.15. For a lower average 
albedo (2%) we find a shallow slope, q — 1.64 ± 0.21, with 
a larger uncertainty as some objects in such a distribution 
will occasionally have Ar ~ 0, giving very large sizes and 
therefore a large spread in slopes. A higher average albedo 
(6%) gives q = 1.89 ± 0.12, a steeper and less uncertain 
slope. A narrower distribution (Ar = 0.04 ± 0.005) gives 
q = 1.93±0.09, while a broad distribution (Ar = 0.04±0.02) 
gives q= 1.72 ±0.19. 



4.4 Correction for shape 

4-4-1 Relating true effective radii and snap-shot radii 

Comets are generally found to be non-spherical, and the 
varying cross-sectional area of a rotating elongated body re- 
flects a changing amount of light, which would imply a differ- 
ent effective radius if a 'snap-shot' measurement was made 
at different rotational phase. Nuclei are often described as 
tri-axial ellipsoids, with semi-major axes a ^ b = c, as this 
is the simplest non-spherical body. Except for comets that 
have been visited by spacecraft, there generally is not suf- 
ficient data to make a more complex model of the shape. 
When producing a size distribution of these bodies, we are 
interested in the effective radius of a sphere with the same 
volume, given by tn = \/ ab 2 , but what we measure us- 
ing equation [2] is the radius corresponding to the projected 
cross sectional area at the time of observation. The pro- 
jected area of a tri-axial ellipsoid, at rotational phase <j> and 
with the pole orien tated at e to the line of sight, is given by 
jLamv et alj|2004 ): 



S(cf>, e) — nab 



IP 



sin 2 e + 



cos 2 e 



IP 



(3) 



The radius found from a snap-shot observation r BB is the 
one corresponding to this area. Substituting S(4>, e) = nr 2 s 
and ab 2 = into equation [3] we find that the true effective 
radius is related to the snap-shot radius by: 



(!) 



a2 
b 



+ cos 



■ 2 2 

sin e + cos e 



(4) 



This equation can be used to find the average relation be- 
tween rN and r ss by integrating over 0^0^ 180° and 
^ e ^ 90°, which gives a correspondence that is a weak 
function of a/b: 



rN = r ss - 



V2 



(5) 



For an a/b — 1.5 this gives tn = 0.97r ss , and for a more ex- 
treme a/b = 3 we find tn = 0.89r ss . This implies that snap- 
shots will, on average, give a very reasonable measurement 
of the effective radius, with the true radius being between 
90-100% of the reported radius from a single snap-shot. 

For individual measurements this does not hold; some- 
times an observation will happen to fall at light curve min- 
imum, and sometimes it will fall at a maximum. Instead of 
applying a statistical correction to the whole database of 



1 1.5 2 2.5 3 

a/b 

Figure 5. Distribution of observed minimum axial ratios for 
cometary nuclei, from the magnitude ranges listed in table [2] 



snap-shots under the assumption that they will average out, 
we use MC simulations to generate the 'true' radius of the 
nucleus from a snap-shot observation using equation [4] We 
do this by selecting a random shape (a/b) and orientation 
(4>, e) for the nucleus, together with the radius r SB from the 
snap-shot magnitude (and all other assumed parameters) 
given by equation [2] Of course the 'true' radius generated 
this way may not match reality for any individual case, but 
over the N MC runs this will generate a distribution of the 
possible effective radii that match the snap-shot photome- 
try. This method has the advantage that it doesn't simply 
apply a correction for snap-shots that shift the radii up or 
down uniformly, but produces a family of CSDs describing 
how the overall distribution is affected by the possible vari- 
ation in individual comet radii due to shape, allowing us to 
measure the uncertainty introduced. 

The MC method requires that we choose appropriate 
distributions for the random parameters. For <j> and e this is 
trivial; any value in the ranges (f> ^ 180° and ^ e ^ 90° 
respectively is equally likely as the rotational phase and pole 
orientation is entirely random for comets, as torques from 
jets can change the spin state in any direction on a short 
timescale. When a light curve is measured, one obtains the 
magnitude averaged over rotational phase, but e is generally 
unknown. In this case we follow the same technique, but set 
= and only randomise e. The choice of axial ratio is 
more difficult, as we need to understand the distribution of 
nucleus shapes. In the next subsection we discuss the current 
knowledge of this distribution. 

4-4.2 Nucleus shapes 

Some constraint can be placed on the shape of an individual 
nucleus from its light curve, as a minimum value of a/b can 
be found from the amplitude of the light curve Am: 

/ a\ nab _ Area ma x Flux max 

V b J min nb 2 Area mi 



Flux n 



= 10 u 



(6) 



The distribution of minimum axial ratios from the light 
curves available in our catalogue is shown in figure it is 
peaked at a/b ~ 1.5 and can be reasonably well described as 
having a decrease in N as a/b increases, with a lack of bod- 
ies with a/b « 1. However, this describes only the observed 
distribution of minimum axial ratios, not the distribution of 
actual shapes. 
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Figure 6. Distributions of apparent minimum a/b for 10,000 'ob- 
servations' of nuclei with fixed axial-ratio at random e. Input val- 
ues of true a/b = 1.2, 1.5, 1.8, 2.2. 



In order to describe how the observed minimum a/b is 
related to the true axial ratio of the nucleus, we substitute 
equation [3] into the ratio of areas given in equation [6] 



S(<t> = 90°) 



■ 2 2 

sin e + cos e 



(7) 



This returns the true a/b when the nucleus is viewed from 
the equator (e = 90°) and no light curve at all (observed 
a/b — 1) from the pole (e = 0°). The form of this equa- 
tion means that the observed minimum a/b increases slowly 
from 1 as e increases from 0, rising steeply towards the true 
value at e > 60° for highly elongated nuclei (see fig. 1 of 
lLamv et al.l (12004 )1. A histogram of the observed minimum 
a/b for a sequence of light curves of an object, taken at ran- 
dom values of e, is bimodal. It has a larger peak at a/b ~ 1 
and a smaller one at the true a/b (fig. [6]). We show his- 
tograms for four different nuclei with axial ratios of 1.2, 1.5, 
1.8 and 2.2. As the true a/b increases, the difference in the 
relative size of the peaks in the histogram also increases; 
for nuclei with elongations > 1.5 a constant ~ 20% of the 
'observed' light curves give a minimum a/b ~ 1. The num- 
ber of observations that return the true axial ratio decreases 
rapidly as the nucleus becomes more elongated, as the prob- 
ability of seeing a more elongated body near enough to end 
on to get an accurate measurement of the true minimum 
area decreases. 

For a collection of nuclei with varying shapes, the distri- 
bution of observed minimum a/b from light curves at ran- 
dom e is a combination of these histograms. We modelled 
this by generating an input distribution of true shapes and 
finding the resulting minimum a/b distribution 'observed' 
in N — 10, 000 light curves. We find that for any input 
distribution of true axial-ratios, the randomisation of pole 
orientation means that light curves return a distribution of 
observed minimum a/b that is strongly peaked at ~ 1 and 
falls off exponentially. Testing a number of input distribu- 
tions (fiat, Gaussian, exponentially decreasing and power 
laws with decreasing and even increasing slope) all pro- 
duced qualitatively similar output distributions, with only 
slight differences in the tail of the distributions visible in 
fig. [3 The best fit power laws (dN/d(a/b) tx (a/b)~ x ) to 



Figure 7. Left: Input 'true' a/b distributions. Right: Output 'ob- 
served' distributions based on random e. From the top: Flat dis- 
tribution between 1 ^ a/b ^ 3; truncated Gaussian distribution 
(< a/b >= 1.5; <r = 0.6); exponential distribution (r = 0.5); 
power law (x = 5.6); power law (x = —5.6). 



the output distributions have x as 7, 7, 11, 15, and 5 respec- 
tively. Therefore it is impossible to constrain the underlying 
shape distribution from the distribution of the minimum 
axial-ratios observed in light curves. 

In reality, there is an observational bias as nuclei with 
a/bmin ~ 1 produce very small amplitude brightness varia- 
tions, making the light curve un-measureable; it is impossi- 
ble to know that the full range of the light curve has been 
observed if peaks and troughs are not recognisable, unless 
the rotation period of the nucleus is independently known. 
The difference between the appearance of the theoretical 
'observed' distributions and that of fig. [5] is explained by 
this lack of reported light curves with Am ~ 0. 
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q=1 .92 



Figure 8. Shape distributi on for asteroids. Le ft: True a/b distri- 
bution from shape models l|Torppa et alj|200Sh . Right: Observed 
a/b distribution from LCDB. 



4-4-3 Assumed distribution and MC results 

As we cannot constrain the true distribution of axial ra- 
tios from the observed distribution, we are forced to assume 
a distribution for our M-C simulations: We choose to take 
asteroid shapes as a reasonable basis. The observed distri- 
bution of minimum axial ratios from 3666 lig ht curves from 
the A steroid Tight Curve Data Base (LCDB: IWarner et al.l 
l2009h shows the same expected shape: A decrease in N as 
a/b increases but with a lack of objects at a/b ~ 1 where 
the light curve cannot be measured (fig. [8|. However, for 
asteroids we can compare this with the distribution of 'true' 
shapes as a significa nt number have detailed shape models. 
iTorppa et all (j200&t ) published a list of simplified descrip- 
tions of 87 models; the distribution of axial ratios from this 
work is best described by a power law with x = 5.6 (fig. 
[8| , although we note that this sample is dominated by large 
asteroids. 

Using this power law (x = 5.6, normalised over the 
range 1 ^ a/b ^ 3) gives q = 2.01 ± 0.07. Changing the 
normalisation to 1 ^ a/b ^ 5 makes no difference. There is 
also only a very small difference (Aq < 0.05) for a range of 
power law distributions of the axial ratio from x = 2.6 to 
8.6, although the uncertainty on the fit increases with de- 
creasing x. This is most likely due to the smaller variation 
in shapes (and hence smaller variation in corrected size) at 
larger x, where a larger proportion of bodies are approxi- 
mately spherical. 

We also tested other shape distributions. Using a Gaus- 
sian distribution of a/b (average 1.5, a = 0.6; a reasonable 
approximation to the observed distribution of minimum a/b 
values) instead gives q — 1.95 ±0.10. Changing the width 
of the distribution has no appreciable effect on q, while in- 
creasing the mean to a/b — 2.0 slightly decreases the slope 
to q — 1.91 ± 0.11. A flat distribution of a/b values between 
1 and 3 gives q — 1.92 ± 0.12, extending it to a/b — 5 gives 
q = 1.87T0.15. In general, the slope of the CSD decreases as 
the distribution includes more elongated bodies, but due to 
the geometrical effects described above, there is little differ- 
ence for any reasonable shape distribution. While we have 
shown that it is impossible to constrain the underlying shape 
distribution from light curve data alone due to random pole 
orientations, the same effect means that snap-shots tend to 
return approximately the correct effective radius without 
any knowledge of the true shape. Therefore the contribu- 
tion to the uncertainty on the CSD due to the unknown 
shape of nuclei is not very large. 
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Figure 9. Probability map for MC run 39, which gives our final 
result by allowing all parameters to vary. The shading shows the 
average shape of the 10,000 CSDs: darker areas are the bins in 
_R-P(rN > R) space where the majority of CSDs passed through, 
lighter areas show the outlying areas explored at the ends of the 
distributions. A dashed line showing the average slope from this 
run is over-plotted. 



4.5 Including all parameters 

A Monte Carlo run allowing all parameters to vary within 
our preferred distributions for each (Gaussians with the ob- 
served mean and a values for An and j3, and an x = 5.6 
power law a/b distribution) represents our best description 
of what we expect the true size distribution to be, and gives 
an uncertainty including all parameters. The fit to this has 
q — 1.92 ± 0.20, with an average cut off in radius at 1.25 
km, and a best fit shallow slope of 0.18 ± 0.03 at sizes be- 
low this. The same run but fixing /3 and An at the 'known' 
values for the subset of comets with measurements gives 
q = 2.01 ± 0.20. Note the relatively large uncertainty when 
including all parameters at once; we take this as a fair as- 
sessment of the true uncertainty on the CSD slope including 
all unknown factors. 



5 DISCUSSION 

Table [S] shows input parameters and results for all the MC 
runs, including all of the ones varying only a single parame- 
ter at a time described in the previous section and also runs 
allowing all parameters to vary. This shows that, relative to 
our reference CSD with q = 1.97, the largest change in slope 
is found when allowing the albedo to vary. This is mostly due 
to the very shallow slope that results from any albedo dis- 
tribution that includes exceptionally dark (Ar ~ 0) nuclei, 
which are unlikely to be realistic. In any case, this demon- 
strates the importance of better constraining the albedo 
distribution of comets, a n important result that the SEPP- 
CoN survey will provide (1 Fernandez et a l. 2008; L owrv et al.l 
120101) . The only other MC runs with a large difference from 
the reference (Aq ^ 0.1) are the extreme shape distribu- 
tions, which are also unlikely to represent reality. We can 
conclude from this that the uncertainty on the input pa- 
rameter distributions actually has only a small effect, as the 
variation in slopes always falls within the typical uncertainty 
found when varying a single parameter at a time. 

Allowing all parameters to vary within the best-fit dis- 
tributions gives a final value of q = 1.92 ± 0.20. This value 
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Table 5. Inputs and results from Monte Carlo runs. 



# 


Mag^ 


Phase function 


Albedo 


n. J7)| 

Shapc£J 


Fix|_ 


N 


d 

<?c, 


p 


Cut-off 


Fixed assumed values: 
















1 


N 


0.035 


0.04 


Spherical 


N 


1 


0.19 


1.97 


1.26 


2 


N 


0.035 


0.04 


Spherical 


Y 


1 


0.19 


2.00 


1.26 


Adjust photometry: 
















3 


Y 


0.035 


0.04 


Spherical 


N 


1000 


0.18 ±0.01 


1.93 ±0.10 


1.24 ±0.07 


Adjust phase 


/ (All tr\^bb\J 1 L. 
















1 


N 


0.025 


0.04 


Spherical 


N 


1 


0.18 


1.94 


1.17 


5 


N 


0.045 


0.04 


Spherical 


N 


1 


0.19 


1.98 


1.35 


6 


N 


0.065 


0.04 


Spherical 


N 


1 


0.20 


1.99 


1.50 


7 


N 


n okq 4- n m r 

\J.\J\)'J _l_ U.U1U 


0.04 


Spherical 


N 


1000 


0.19 ±0.01 


1.95 ±0.07 


1.38 ±0.07 


8 


N 


0.053 ± 0.005 


0.04 


Spherical 


K 


1000 


0.20 ±0.00 


1.99 ±0.02 


1.42 ±0.03 


9 


N 


n nKQ _i_ n no"; 

W.UUtJ _l_ U.UinJ 


0.04 


Spherical 


N 


1000 


0.19 ±0.02 


1.91 ±0.11 


1.35 ±0.09 


10 


N 


n r\t\o J. nnifj 
u.uuo r_n u.uiu 


0.04 


Spherical 


Y 


1000 


0.19 ±0.01 


2.05 ±0.06 


1.39 ±0.06 


11 


N 


010 4- D 01 

U.U-LU _l_ U.UIU 


0.04 


Spherical 


N 


1000 


0.15 ±0.01 


1.72 ±0.06 


1.02 ±0.06 


12 


N 


090 -1- 01 

U.U-^U □! U.UIU 


0.04 


Spherical 


N 


1000 


0.16 ±0.02 


1.81 ±0.11 


1.13 ±0.11 


13 


N 


0^0 4- 01 

yj.yjrjyj _i_ u.uiu 


0.04 


Spherical 


N 


1000 


0.18 ±0.01 


1.94 ±0.09 


1.26 ±0.09 


14 


N 


040 -1- 01 

U.UIU □! U.UIU 


0.04 


Spherical 


N 


1000 


0.19 ±0.01 


1.97 ±0.05 


1.32 ±0.06 


15 


N 


n^Q + o 01 

U.UOU _1_ u.ulU 


0.04 


Spherical 


N 


1000 


0.19 ±0.01 


1.97 ±0.04 


1.38 ±0.06 


16 


N 


060 + 010 


0.04 


Spherical 


N 


1000 


0.20 ±0.01 


1.97 ±0.05 


1.43 ±0.06 


17 


N 


n (i7n + n oi o 


0.04 


Spherical 


N 


1000 


0.20 ±0.01 


1.95 ±0.05 


1.47 ±0.06 


18 


N 


080 ± 010 


0.04 


Spherical 


N 


1000 


0.20 ±0.01 


1.92 ±0.05 


1.49 ±0.05 


19 


N 


n nqn + n 01 

U.UC7U _1_ U.UIU 


0.04 


Spherical 


N 


1000 


0.20 ±0.01 


1.90 ±0.04 


1.54 ±0.04 


Adjust albedo 


















20 


N 


0.035 


0.044 ± 0.013 


Spherical 


N 


1000 


0.17 ±0.02 


1.85 ±0.15 


1.16 ±0.10 


21 


N 


0.035 


0.020 ± 0.013 


Spherical 


N 


1000 


0.16 ±0.03 


1.64 ±0.21 


1.64 ±0.21 


22 


N 


0.035 


0.060 ± 0.013 


Spherical 


N 


1000 


0.18 ±0.02 


1.89 ±0.12 


1.01 ±0.08 


23 


N 


0.035 


0.040 ± 0.005 


Spherical 


N 


1000 


0.18 ±0.01 


1.93 ±0.09 


1.26 ±0.08 


24 


N 


0.035 


0.040 ± 0.020 


Spherical 


N 


1000 


0.16 ±0.03 


1.72 ±0.19 


1.18 ±0.13 


25 


N 


0.035 


0.044 ± 0.013 


Spherical 


Y 


1000 


0.17 ±0.02 


1.85 ±0.14 


1.16 ±0.10 


Adjust shape: 


















26 


N 


0.035 


0.04 


Flat, 1 ^ a/b ^ 3 


N 


1000 


0.18 ±0.02 


1.92 ±0.12 


1.17 ±0.09 


27 


N 


0.035 


0.04 


Flat, 1 ^ a/b ^ 5 


N 


1000 


0.17 ±0.02 


1.87 ±0.15 


1.09 ±0.09 


28 


N 


0.035 


0.04 


Gaussian, 1.5 ± 0.6 


N 


1000 


0.18 ±0.01 


1.95 ±0.10 


1.22 ±0.08 


20 


N 


0.035 


0.04 


Gaussian, 1.5 ± 0.3 


N 


1000 


0.18 ±0.01 


1.96 ±0.08 


1.23 ±0.07 


30 


N 


0.035 


0.04 


Gaussian, 1.5 ± 0.9 


N 


1000 


0.18 ±0.01 


1.94 ±0.11 


1.20 ±0.08 


31 


N 


0.035 


0.04 


Gaussian, 2 ± 0.6 


N 


1000 


0.18 ±0.01 


1.91 ±0.11 


1.18 ±0.08 


32 


N 


0.035 


0.04 


Power law, x = -5.6, 


N 


1000 


0.18 ±0.01 


2.01 ±0.07 


1.26 ±0.06 










1 ^ a/b ^ 3 












33 


N 


0.035 


0.04 


Power law, x = -5.6, 


N 


1000 


0.18 ±0.01 


2.02 ±0.07 


1.26 ±0.07 










U a/6 sj 5 












34 


N 


0.035 


0.04 


Power law, x = -4.6, 


N 


1000 


0.18 ±0.01 


2.00 ±0.08 


1.25 ±0.07 










Ha/K5 












35 


N 


0.035 


0.04 


Power law, x = -6.6, 


N 


1000 


0.18 ±0.01 


2.02 ±0.06 


1.27 ±0.06 










1 ^ a/b < 5 












36 


N 


0.035 


0.04 


Power law, x = -2.6, 


N 


1000 


0.18 ±0.01 


1.97 ±0.11 


1.22 ±0.08 










U «/K 5 












37 


N 


0.035 


0.04 


Power law, x = -8.6, 


N 


1000 


0.18 ±0.00 


2.02 ±0.04 


1.27 ±0.05 










U a/6 sj 5 












Adjust all parameters: 
















38 


Y 


0.035 ± 0.020 


0.040 ± 0.040 


Gaussian, 1.5 ± 0.6 


N 


10000 


0.16 ±0.04 


1.51 ±0.22 


1.17 ±0.19 


39 


Y 


0.053 ± 0.016 


0.044 ± 0.013 


Power law, x = -5.6, 


N 


10000 


0.18 ±0.03 


1.92 ±0.20 


1.25 ±0.13 










1 ^ a/b sC 3 












10 


Y 


0.053 ± 0.016 


0.044 ± 0.013 


Power law, x = -5.6, 


Y 


10000 


0.19 ±0.02 


2.01 ±0.20 


1.28 ±0.12 



1 < a/K 3 



a Vary input magnitude within photometric error bars? 
^ Assumed shape distribution for nuclei. 

c Use fixed values for /3 and An when known for a given comet? If no, all comets have values selected from the specified 
distribution, even those where these values are independently constrained. 

" qi is the best fit at r'N ^ cut-off. 

e 52 is the best fit at tn > cut-off, i.e. the quoted q. 



12 C. Snodgrass et al. 



is remarkably close to our reference value, again implying 
that the combined effect of the various uncertain parame- 
ters is actually relatively small. This helps to validate the 
assumptions made in previous CSD estimates, but the real 
significance of this work is that the uncertainty on this value 
has been rigourously determined by including all individual 
uncertainties on each radius. Figure [9] shows how the CSDs 
varied in this MC simulation. 

Comparing this with the previous CSD estimates listed 
in Table [T] shows that the value we find for q is in a similar 
range to recent results, if slightly steeper than the consensus 
value, and when we consider the true uncertainty we find the 
majority of results are encompassed within 2a. While pre- 
vious works assumed constant values for An and /3 and did 
not consider the uncertainties on these, we show that the 
CSD slope including all sources of uncertainty agrees with 
these earlier results. It is worth noting that our testing of 
various distributions for the assumed parameters never pro- 
duced a slope as steep as th e considerably l a rger q found by 
iFernandez et al.l \l99$ ) and lTancredi et al.l (|2006l) . even for 
extreme distributions we regard as unlikely to be realistic. 
We suspect that the difference is due to the higher accuracy 
of the photometric measurements used in our catalogue, all 
of which come from large aperture telescopes. It is possible 
that a similar MC analysis of th e uncertainty on q due to the 
photometric uncertainties in the lTancredi et all (|2006h cata- 
logue would give a large error bar, making it consistent with 
the other results, however this catalogue does not include all 
the necessary information (in its currently published form) 
to apply our technique. 

The CSD slope we find is very close to the theoreti - 
cal value of q = 2.04 found bv lO'Brien fe Greenberel (|2003T l 
for a collisionally relaxed population of strengthless bodies. 
However, the large uncertainty we find means that we can- 
not claim that this is evidence that JFC nuclei came from 
collisional disruption of larger bodies; we can simply state 
that the CSD slope is consistent with that interpretation. It 
is important to remember that the activity of comets alters 
their sizes as they lose significant amounts of material with 
each perihelion passage, and therefore we would not neces- 
sarily expect to see a CSD that matched a collisional one, 
even if the comets are collisional fragments. Both 'normal' 
cometary activity and major mass loss events (nucleus split- 
ting or outbursts) are highly variable between comets and 
from one orbit to another, so the change in size due to mass 
loss from a nucleus is a highly complicated process. More 
detailed theoretical models describing the expected CSD of 
nuclei for various formation (and evolution) scenarios are re- 
quired for comparison with observations, however the large 
uncertainty on q found by this paper shows that we must 
first provide better observational constraints on all physi- 
cal properties of nuclei (albedos, phase functions and shape 
models) to produce observed CSDs with sufficient accuracy 
to compare with these models. 
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